package MathLib.MyMatrices;

public class CEignValues extends CMatrix {
	protected static int maxIterNum = 1000000;
    protected double eps;
//    	protected double[][] R;
//    	protected double[][] Q;
//		protected double[][] L;
//		protected double[][] U;
    
    
    public CEignValues(){
    	super();

    	this.eps = 0.00000001;
    } 
    
    public CEignValues(int n,double eps){
    	super(n);

    	this.eps = eps;
    }
   
    public CEignValues(CMatrix A,double eps) throws MatrixEception{
    	super(A);
    	if (A.n!=A.m) throw new MatrixEception(3);    	
    	this.eps = eps;
    }
    
    public CEignValues(CMatrix A) throws MatrixEception{
    	super(A);
    	if (A.n!=A.m) throw new MatrixEception(3);    	
    	this.eps = 0.00000001;
    }
    
    public double PMMethod(double[] y,double delta) throws MatrixEception {
    	if(y.length!=n){
    		throw new MatrixEception(7);
    	}
    	
    	double norma = 0.0;
    	for(int i = 0; i<n;i++){
    		norma+=y[i]*y[i];
    	}
    	norma = Math.sqrt(norma);
    	double[] x = new double[n];
    	
    	for(int i = 0; i<n;i++){
    		x[i] = y[i]/norma;
    	}
		int k = 0;
		double l1=0.0,l2;
    	while(k<maxIterNum){
    		k++;
    		//2
    		for(int i = 0; i<n;i++){
    			y[i]=0.0;
    			for(int j = 0; j<n;j++){
    	    		y[i]+= x[j]*Matrix[i][j];
    	    	}
        	}
    		//3
    		for(int i = 0; i<n;i++){
        		norma+=y[i]*y[i];
        	}
        	norma = Math.sqrt(norma);
        	
        	for(int i = 0; i<n;i++){
        		x[i] = y[i]/norma;
        	}
    		//4
        	l2=l1;
        	l1=0.0;
        	int p=0;
        	for(int i = 0; i<n;i++){
        		if(Math.abs(x[i])>delta){
        			l1+=y[i]/x[i];
        			p++;
        		}
        	}
        	l1/=p;
        	if((k!=0)&&(Math.abs(l1-l2)<eps)){
        		y=x;
        		return l1;        		
        	}
    		
    	}
    	y=x;
    	return l1;
	}
    
    public double SPMethod(double[] y) throws MatrixEception {
    	if(y.length!=n){
    		throw new MatrixEception(7);
    	}
    	
    	double norma = 0.0;
    	for(int i = 0; i<n;i++){
    		norma+=y[i]*y[i];
    	}
    	norma = Math.sqrt(norma);
    	double[] x = new double[n];
    	
    	for(int i = 0; i<n;i++){
    		x[i] = y[i]/norma;
    	}
    	
		int k = 0;
		double l1=0.0,l2;
    	while(k<maxIterNum){
    		k++;
    		//2
    		for(int i = 0; i<n;i++){
    			y[i]=0.0;
    			for(int j = 0; j<n;j++){
    	    		y[i]+= x[j]*Matrix[i][j];
    	    	}
        	}
    		//3
    		double norma2 = 0.0;
    		norma = 0.0;
    		for(int i = 0; i<n;i++){
        		norma+=y[i]*y[i];
        		norma2+=y[i]*x[i];
        	}
    		norma = Math.sqrt(norma);
    		norma2 = Math.sqrt(norma2);
    		l2=l1;
    		l1=norma/norma2;
        	
        	for(int i = 0; i<n;i++){
        		x[i] = y[i]/norma;
        	}
    		//4

        	if((Math.abs(l1-l2)<eps)){
        		return l1;        		
        	}
    		
    	}
    	return l1;
	}

    public void JakobiRotation(double[] values) {
        int num = n*n-n;
        int[] imas = new int[2];
    	for (int k = 1;k < num;k++){
            int i=0,j=0;
            double TmpMax = vedychuj_element(imas);
            i = imas[0];
            j = imas[1];
            //if(Math.abs(TmpMax)<eps){break;}
            double p = 2 * TmpMax;
            double q = Matrix[i][i] - Matrix[j][j];
            double d = Math.sqrt(p*p + q*q);

            double r,c,s;
            if (Math.abs(q)>0.000000001){
                r = Math.abs(q)/(2*d);
                c = Math.sqrt(0.5 + r);
                s = Math.sqrt(0.5 - r)*Math.signum(p*q);
            }
            else{
                 c = Math.sqrt(2)/2;
                 s = c;
            }

         // diagonal elements
         double temp1 = Matrix[i][i],temp2 = Matrix[j][j];
         Matrix[i][i] = c*c*temp1 + s*s*temp2 + 2*c*s*TmpMax ;
         Matrix[j][j] = s*s*temp1 + c*c*temp2 - 2*c*s*TmpMax ;
         // set to Zero maximal element
         Matrix[i][j] = 0;
         Matrix[j][i] = 0;

         for( int i1 = 0; i1<n;i1++){
            if ((i1!=j) && (i1!=i)) {
                Matrix[i][i1] = c*Matrix[i1][i] + s*Matrix[i1][j];
                Matrix[i1][i] = Matrix[i][i1];
                Matrix[j][i1] = -s*Matrix[i1][i] + c*Matrix[i1][j];
                Matrix[i1][j] = Matrix[j][i1];
            }
            }
        }

        for( int i1 = 0; i1<n;i1++){
            values[i1] = Matrix[i1][i1];
        }


	}
    
    private double vedychuj_element(int[] imas){
    	int i1=n-1,j1=0;
    	double temp = Math.abs(Matrix[i1][j1]);
    	for(int i=0;i<n;i++){
    		for(int j=i+1;j<n;j++){
    			if((i!=j)&&(Math.abs(Matrix[i][j])>temp)){
    				temp = Math.abs(Matrix[i][j]);
    				i1=i;
    				j1=j;
    			}
        	}
    	}
        imas[0] = i1;
        imas[1] = j1;
        return temp;
    }
}
